rm(list=ls())


library(lattice)
library(gridExtra)
setwd("~/Dropbox/Rule of Law (EU Enlargement)")

source("compliance function July 23 2012.R")

# 100 Rounds; N=40; threshold = 0.75
nrounds <- 100
n<-40
t<-0.75

# True population support is high --- 90% (e.g. 36 1s and 4 0s)
pop <- c(rep(1,times=36),rep(0,times=4))

probs <- seq(0,1,by=0.02)
parmat<-expand.grid(probs,probs)
parmat<-cbind(parmat,rep(NA,nrow(parmat)))
colnames(parmat)<-c("pvec","rvec","prop.out")


for (j in 1:nrow(parmat)){

p <- parmat[j,1]
r <- parmat[j,2]
q <- 0.5

post <- rep(NA,nrounds+1)
post[1] <- p

for (i in 1:nrounds){
	a <- comp(p,n,pop,t,r,q)
	p <- mean(a$out)
	post[i+1] <- p	
	meanpac <- mean(a$pac)
}

parmat[j,3] <- meanpac

}

write.csv(parmat,file="~/Dropbox/Rule of Law (EU Enlargement)/parmatpop90.csv")


#########################################################


rm(list=ls())

setwd("~/Dropbox/Rule of Law (EU Enlargement)")

source("compliance function July 23 2012.R")

# 100 Rounds; N=40; threshold = 0.75
nrounds <- 100
n<-40
t<-0.75

# True population support is low --- 10% (e.g. 4 1s and 36 0s)
pop <- c(rep(1,times=4),rep(0,times=36))

probs <- seq(0,1,by=0.02)
parmat<-expand.grid(probs,probs)
parmat<-cbind(parmat,rep(NA,nrow(parmat)))
colnames(parmat)<-c("pvec","qvec","prop.out")

for (j in 1: nrow(parmat)){

p <- parmat[j,1]
q <- parmat[j,2]
r <- 0.5

post <- rep(NA,nrounds+1)
post[1] <- p

for (i in 1:nrounds){
	a <- comp(p,n,pop,t,r,q)
	p <- mean(a$out)
	post[i+1] <- p	
	meanpac <- mean(a$pac)
}

parmat[j,3] <- meanpac

}

write.csv(parmat,file="~/Dropbox/Rule of Law (EU Enlargement)/parmatpop10.csv")


#######################################################

rm(list=ls())

setwd("~/Dropbox/Rule of Law (EU Enlargement)")

source("compliance function July 23 2012.R")

# 100 Rounds; N=40; threshold = 0.75
nrounds <- 100
n<-40
t<-0.75

# True population support is middle --- 50% (e.g. 20 1s and 20 0s)
pop <- c(rep(1,times=20),rep(0,times=20))

probs <- seq(0,1,by=0.03)
parmat<-expand.grid(probs,probs,probs)
parmat<-cbind(parmat,rep(NA,nrow(parmat)))
colnames(parmat)<-c("pvec","rvec","qvec","prop.out")

for (j in 1: nrow(parmat)){

p <- parmat[j,1]
r <- parmat[j,2]
q <- parmat[j,3]

post <- rep(NA,nrounds+1)
post[1] <- p

for (i in 1:nrounds){
	a <- comp(p,n,pop,t,r,q)
	p <- mean(a$out)
	post[i+1] <- p	
	meanpac <- mean(a$pac)
}

parmat[j,4] <- meanpac

}

write.csv(parmat,file="~/Dropbox/Rule of Law (EU Enlargement)/parmatpop50.csv")



######## Plots ######

parmat90<-read.csv("~/Dropbox/Rule of Law (EU Enlargement)/parmatpop90.csv")
parmat10<-read.csv("~/Dropbox/Rule of Law (EU Enlargement)/parmatpop10.csv")
parmat50<-read.csv("~/Dropbox/Rule of Law (EU Enlargement)/parmatpop50.csv")

trellis.par.set("axis.line",list(col="transparent"))

# Plots for pop=90
pdf(file="~/Dropbox/Rule of Law (EU Enlargement)/pop90.pdf",width=15,height=10)
h1= histogram(parmat90$prop.out,
			col="gray",
			ylab=list(label="Percent Simulations"),
			xlab=list(label="Percent Population Obeying Law"),   
			)
wf1= wireframe(parmat90$prop.out~parmat90$pvec+parmat90$rvec,
	xlab="p",
	ylab="r",
	zlab=list(label="Proportion Obeying",rot=90),
	drape=T,
	colorkey=F,
	scales =list(arrows=F,cex=.45, col="black", font=3, tck=1))
grid.arrange(h1,wf1,ncol=2)
dev.off()

# Plots for pop=10
pdf(file="~/Dropbox/Rule of Law (EU Enlargement)/pop10.pdf",width=15,height=10)
h2= histogram(parmat10$prop.out,
			col="gray",
			ylab=list(label="Percent Simulations"),
			xlab=list(label="Percent Population Obeying Law"),   
			)
wf2= wireframe(parmat10$prop.out~parmat10$pvec+parmat10$qvec,
	xlab="p",
	ylab="q",
	zlab=list(label="Proportion Obeying",rot=90),
	drape=T,
	colorkey=F,
	scales =list(arrows=F,cex=.45, col="black", font=3, tck=1))

grid.arrange(h2,wf2,ncol=2)
dev.off()


# Plots for pop=50
pdf(file="~/Dropbox/Rule of Law (EU Enlargement)/pop50.pdf",width=12,height=10)
h3= histogram(parmat50$prop.out,
			col="gray",
			ylab=list(label="Percent Simulations"),
			xlab=list(label="Percent Population Obeying Law"),   
			)
wf3= wireframe(parmat50$prop.out~parmat50$pvec+parmat50$rvec,
	xlab="p",
	ylab="r",
	zlab=list(label="Proportion Obeying",rot=90),
	drape=T,
	colorkey=F,
	scales =list(arrows=F,cex=.45, col="black", font=3, tck=1))

wf4= wireframe(parmat50$prop.out~parmat50$pvec+parmat50$qvec,
	xlab="p",
	ylab="q",
	zlab=list(label="Proportion Obeying",rot=90),
	drape=T,
	colorkey=F,
	scales =list(arrows=F,cex=.45, col="black", font=3, tck=1))

wf5= wireframe(parmat50$prop.out~parmat50$rvec+parmat50$qvec,
	xlab="r",
	ylab="q",
	zlab=list(label="Proportion Obeying",rot=90),
	drape=T,
	colorkey=F,
	scales =list(arrows=F,cex=.45, col="black", font=3, tck=1))

grid.arrange(h3,wf3,wf4,wf5,ncol=2)
dev.off()


wireframe(parmat$prop.out~parmat$pvec+parmat$qvec)
wireframe(parmat$prop.out~parmat$pvec+parmat$qvec)
wireframe(parmat$prop.out~parmat$pvec+parmat$rvec)
